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Abstract 

I discuss liow to sample the space of parton momenta in order to best perform 
the numerical integrations that lead to a calculation of three jet cross sections 
and similar observables in electron-positron annihilation. 
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I. INTRODUCTION 



There is an important class of computer programs that do calculations in quantum 
chromodynamics (QCD) in which the calculation is performed at next-to-leading order in 
perturbation theory and allows for the determination of a variety of characteristics of the 
final state. This paper concerns a significant technical issue that arises in such programs 
when a "completely numerical" integration algorithm is used: how should one choose the 
integration points? 

I consider the calculation of "three-jet-like" observables in e'^e" annihilation. A given 
program in the class considered can be used to calculate a jet cross section (with any infrared 
safe choice of jet definition) or observables like the thrust distribution. Such a program 
generates random partonic events consisting of three or four final state quarks, antiquarks, 
and gluons. Each event comes with a calculated weight. A separate routine then calculates 
the contribution to the desired observable for each event, averaging over the events with 
their weights. 

The weights are treated as probabilities. However, these weights can be both positive 
or negative. This is an almost inevitable consequence of quantum mechanics. The calcu- 
lated observable is proportional to the square of a quantum amplitude and is thus positive. 
However, as soon as one divides the amplitude into pieces for purposes of calculation, one 
finds that while the square of each piece is positive, the interference terms between different 
pieces can have either sign. Thus the kind of program discussed here stands in contrast 
to the tree-level event generators in which, by simplifying the physics, one can generally 
arrange to have all the weights be positive, or, even, to be all equal to 1. 

To understand the algorithms used in the class of programs described above, it is best to 
think of the programs as performing integrations over momenta in which the quantum matrix 
elements and the measurement functions form the integrand. There are two basic algorithms 
for performing the integrations. The older is due to Ellis, Ross, and Terrano (ERT) (The 
first general purpose implementation of this method for three-jet-like observables in e^e~ 
annihilation was that of Kunszt and Nason 0). In this method, some of the integrations 
are performed analytically ahead of time. The other integrations are performed numerically 
by the Monte Carlo method. The integrations are divergent and are regulated by analytical 
continuation to 3 — 2e space dimensions and a scheme of subtractions or cutoffs. The second 
method is much newer In this method, all of the integrations are done by Monte Carlo 
numerical integration. With this method, the integrals are all convergent (after removal of 
the ultraviolet divergences by a straightforward renormalization procedure). 

Since this method is quite new, one cannot yet say for what problems it might do better 
than the now standard ERT method. I can point out that the current incarnation of the 
numerical method has convergence problems for three jet quantities that receive important 
contributions from final states that are close to the two jet limit. For example, one gets 
good results for the three jet cross section or for the thrust distribution away from T = 1, 
but poorly converging results for the energy-energy correlation function. On the other hand, 
the numerical method offers evident advantages in flexibility to modify the integrand. Thus 
one should be able to attack problems that are not accessible to the ERT method. 

The numerical integration method exists as computer code with accompanying tech- 
nical notes and many of the basic ideas behind it have been described in the two papers 
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The present paper is devoted to an exposition of the choice of integration points needed 
for this method. This subject may seem unimportant. Furthermore, it may seem to be a just 
a part of the numerical algorithms branch of computer science and thus to be uninteresting 
from the point of view of physics. However, practitioners of the art of numerical integrations 
in physics know that the choice of integration points is extremely important. If everything 
else is perfect but the integration points are badly chosen, the computer will, like Sisyphus, 
be faced with an unending task because the precision of the answer will not improve as the 
computational time increases. Furthermore, the problem is a problem in physics. All that 
the theory of numerical integration tells us is that the density of integration points should be 
well matched to the structure of the function to be integrated. The function to be integrated 
comes from Feynman diagrams, so we are really faced with the problem of understanding the 
structure of the quantum scattering processes. Perhaps surprisingly, the aspects of quantum 
scattering that are important in numerical method are completely different than the aspects 
that are important in the ERT method. 

In Sec. II below, I briefly review the basics of the numerical integration method with 
the aim of setting the notation and making it possible to read this paper independently of 
the previous two papers. Then, in the sections that follow, I turn to the main issue of this 
paper, the choice of integration points. I try to keep this discussion succinct. The main 
features of the method are covered, but implementation details are left to the program code 
and its accompanying technical notes P,^. 

The exposition begins in Sec. Ill with a review of the general principles of Monte Carlo 
integration as they apply to this calculation. Sec. IV is the heart of the paper and covers 
the organization of the method used to choose the most integration points in the most 
important integration regions. Part of the method used considers the possible three parton 
final states. This part is explained in Sec. V. The other parts of the method concern the 
singularities of virtual graphs and their relation to the energy denominators in 2 parton to 
n parton scattering. The various possibilities are covered in Sees. VI through IX. Of special 
importance is the elliptical coordinate system introduced in Sec. VI. A brief summary is 
presented in Sec. X. 



II. REVIEW OF THE NUMERICAL METHOD 

Let us begin with a precise statement of the problem. We consider an infrared safe three- 
jet-like observable in e~^e" hadrons such as a particular moment of the thrust distribution. 
The observable can be expanded in powers of as/vr. 



^aM, aMocK/Trf- (1) 



The order a1 contribution has the form 



dpidp2dp3 



1 f da^"^^ 
+ — / dpidp2dp3 i Ss{pi,p2,P3) (2) 
61 J drndvodn?. 
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FIG. 1. Two cuts of one of the Feynman diagrams that contribute to e^e 



hadrons . 



+ 77 



dpidp2dp^dpi 




dpidp2dpsdp4 



Here the da^^^ are the order contributions to the parton level cross section, calculated with 
zero quark masses. Each contains momentum and energy conserving delta functions. The 
dcrlfl include ultraviolet renormalization in the MS scheme. The functions S describe the 
measurable quantity to be calculated. We wish to calculate a "three-jet-hke" quantity. That 
is, ^2 = 0. The normalization is such that iS„ = 1 for n = 2,3,4 would give the order 
perturbative contribution the total cross section. There are, of course, infrared divergences 
associated with Eq. (^. For now, we may simply suppose that an infrared cutoff has been 
supplied. 

The measurement, as specified by the functions Sn, is to be infrared safe, as described 
in Ref. 0: the iS„ are smooth functions of the parton momenta and 



for < A < 1. That is, collinear splittings and soft particles do not affect the measurement. 

It is convenient to calculate a quantity that is dimensionless. Let the functions Sn be 
dimensionless and eliminate the remaining dimensionality in the problem by dividing by (Tq, 
the total e~^e~ cross section at the Born level. Let us also remove the factor of (as/vr)^. 
Thus, we calculate 



Let us now see how to set up the calculation of X in a convenient form. We note that X is 
a function of the cm. energy ^/s and the MS renormalization scale /i. We will choose /i to be 
proportional to i/i: fi = AuvV^- Then X depends on A. But, because it is dimensionless, 
it is independent of y/s. This allows us to write 



<Sn+l{Pl, >^Pn, (1 - A)p„) = Sn{pi, ■ ■ ■ , Pn 



(3) 



X 



(4) 



(To (as/ir)'^' 




(5) 



where h is any function with 




(6) 
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The quantity X can be expressed in terms of cut Feynman diagrams, as in Fig. ^ The 
dots where the parton hues cross the cut represent the function Sn{pi, • • • , Pn)- Each diagram 
is a three loop diagram, so we have integrations over loop momenta , I2 and li^. We first 
perform the energy integrations. For the graphs in which four parton lines cross the cut, 
there are four mass-shell delta functions S{kj). These delta functions eliminate the three 
energy integrals over l^, I2, and /g as well as the integral (P) over ^/s. For the graphs in 
which three parton lines cross the cut, we can eliminate the integration over y/s and two of 
the Zj integrals. One integral over the energy E in the virtual loop remains. We perform 
this integration by closing the integration contour in the lower half E plane. This gives a 
sum of terms obtained from the original integrand by some simple algebraic substitutions. 
Having performed the energy integrations, we are left with an integral of the form 

I = Y. f dh dl2 dh E 9{G, C; k, I2, h). (7) 

G C 

Here there is a sum over graphs G (of which one is shown in Fig. |I]) and there is a sum 
over the possible cuts C of a given graph. The problem of calculating X is now set up in a 
convenient form for calculation. 

If we were using the Ellis-Ross- Terrano method, we would put the sum over cuts outside 
of the integrals in Eq. (|^. For those cuts C that have three partons in the final state, 
there is a virtual loop. We can arrange that one of the loop momenta, say /i, goes around 
this virtual loop. The essence of the ERT method is to perform the integration over the 
virtual loop momentum analytically ahead of time. The integration is often ultraviolet 
divergent, but the ultraviolet divergence is easily removed by a renormalization subtraction. 
The integration is also typically infrared divergent. This divergence is regulated by working 
in 3 — 2e space dimensions and then taking e — > while dropping the 1/e" contributions 
(after proving that they cancel against other contributions). After the li integration has 
been performed analytically, the integrations over I2 and Z3 can be performed numerically. 
For the cuts C that have four partons in the final state, there are also infrared divergences. 
One uses either a "phase space slicing" or a "subtraction" procedure to get rid of these 
divergences, cancelling the 1/e" pieces against the 1/e" pieces from the virtual graphs. In 
the end, we are left with an integral / dli dl2 dl^ in exactly three space dimensions that can 
be performed numerically. 

In the numerical method, we keep the sum over cuts C inside the integrations. We take 
care of the ultraviolet divergences by simple renormalization subtractions on the integrand. 
We make certain deformations on the integration contours so as to keep away from poles of 
the form l/[Ep — Ej ± ie], where Ep is the energy of the final state and Ej is the energy 
of an intermediate state. Then the integrals are all convergent and we calculate them by 
Monte Carlo numerical integration. 

Let us now look at the contour deformation in a little more detail. We denote the 
momenta {h,h,h} collectively by I whenever we do not need a more detailed description. 
Thus 

^=Y.[dl T.9{G,C;l). (8) 

G C 
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For cuts C that leave a virtual loop integration, there are singularities in the integrand of 
the form Ep — Ej + ie (or Ep — Ej — ie if the loop is in the complex conjugate amplitude 
to the right of the cut). Here Ep is the energy of the final state defined by the cut C 
and Ej is the energy of a possible intermediate state. We will examine the nature of these 
"scattering singularities" in some detail in the next section. For now, all we need to know 
is that Ep — Ej = on a surface in the space of h for fixed I2 and I3 if we pick li to be the 
momentum that flows around the virtual loop. These singularities do not create divergences. 
The Feynman rules provide us with the ie prescriptions that tell us what to do about the 
singularities: we should deform the integration contour into the complex li space so as to 
keep away from them. Thus we write our integral in the form 



Here in is a purely imaginary nine-dimensional vector that we add to the real nine- 
dimensional vector / to make a complex nine-dimensional vector. The imaginary part k 
depends on the real part /, so that when we integrate over /, the complex vector / + in lies 
on a surface, the integration contour, that is moved away from the real subspace. When 
we thus deform the contour, we supply a jacobian J' = det{d{l + iK)/dl). (See Ref. for 
details.) 

The amount of deformation k, depends on the graph G and, more significantly, the cut C. 
For cuts C that leave no virtual loop, each of the momenta li, I2, and flows through the 
final state. For practical reasons, we want the final state momenta to be real. Thus we set 
K = for cuts C that leave no virtual loop. On the other hand, when the cut C does leave a 
virtual loop, we choose a non-zero k. We must, however, be careful. When k = there are 
singularities in g on certain surfaces that correspond to collinear parton momenta. These 
singularities cancel between g for one cut C and g for another. This cancellation would be 
destroyed if, for / approaching the collinear singularity, k = for one of these cuts but not 
for the other. For this reason, we insist that for all cuts (7, k ^ as / approaches one of 
the collinear singularities. The details can be found in Ref. 0]. All that is important here 
is that K — i> quadratically with the distance to a collinear singularity. 

Much has been left out in this brief overview, but we should now have enough background 
to see the requirements for a sensible choice of integration points. 

III. GENERAL PRINCIPLES FOR THE CHOICE OF POINTS 

In the following sections, we discuss the choice of integration points for the evaluation of a 
given graph. In this section, we summarize the general principles of Monte Carlo integration 
as they apply to this calculation. 

We wish to perform an integral of the form 



^ = Y. h^Y. ^(^' C*; /) g{G, C] I + m(G, C; /)). 



(9) 



G C 




G 



(10) 



where 
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f{G- /) = 3i |E JiG, C; I) 9{G, C; I + m{G, C; I)) 



(11) 



(We can take the real part because we know in advance that 1 is real.) Using the Monte Carlo 
technique, for each graph G we choose a large number i{G)N of points I, with J2g^{G) — 1 
so that the total number of points is A^. We sample the points / for graph G at random with 
a density p{G; I), with normalization / dl p{G; I) = 1. Then the integral is approximated by 



G 



(12) 



This is an approximation for the integral in the sense that, if we repeat the procedure a 
lot of times, the expectation value for is 



The expected r.m.s. error is S, where 

With a little manipulation, one can rewrite this as 



N 



Ee(G) 



G 



A(G) 



EA(^)]'+(EA(i^)) 



(13) 



(14) 



(15) 



where 



A{Gf 



dl p{G; I) 



\fiG-l)\ 



dl' \f{G;l') 



+ 



p{G;l) 

J dl (|/(G;0| + /(G;0) x J dl' {\f{G;l')\-f{G;l')). 



(16) 



We see, first of all, that the expected error decreases proportionally to 1/viV. Second, we 
see that for a given choice of the density functions p(G;/), the ideal choice of the factors 
^(G) is ^(G) oc A(G). This is, in fact, easy to implement. Third, the ideal choice of p{G;l) 
is p{G; I) oc \ f{G; l)\. This is, in fact, impossible to implement. 

Although it is not possible to choose p{G;l) oc |/(G;/)|, at least we can choose it so 
that \f{G;l)\/p{G;l) is not singular at the singularities of |/(/)|. Furthermore, we can try 
to make p{G; I) big at places where we know that |/(G; /)| is big. 

We will build the general samphng method out of elementary samphng methods. That 
is, we will find a number of simple methods to choose points / for our graph. Let the density 
of points for the ith elementary sampling method be Pi{G; I) (normalized to / dl Pi{G; I) = 1). 
Then we devote a fraction Xi{G) of the points to the choice with density Pi{G; I) and obtain 
a net density 



p{G;l)^^K{G)pi{G;l). 



(17) 



7 



In this way, we make the samphng problem manageable. If we know that \f{G; l) \ is big near 
a certain point or surface in the space of loop momenta, we can design one of the elementary 
sampling methods so that the corresponding Pi{G; I) is big there. In undertaking this task, 
we do not have to simultaneously arrange that Pi{G; I) be big at other places where \f{G; l)\ 
is big. 

In the following sections, we discuss the elementary sampling methods. We imagine that 
we are dealing with only one specific graph G, so we suppress the index G in the notation. 

IV. ORGANIZATION OF THE SAMPLING METHOD 

We consider a fixed (uncut) graph, dropping references to the label G of the graph. 
As mentioned in the previous section, we sample points / in the space of loop momenta 
according to several elementary sampling methods, each labelled by and index i and having 
a density of points Pi(l). The net density is then p(/) = J2KPi{l)- 

We first address the identification of loop momenta. There are eight propagator momenta 
kp. (See, for example. Fig. |T].) The three loop momenta /l can in general be any three linearly 
independent linear combinations of the kp. We will keep the choice simple by choosing three 
of the kp to be the loop momenta. However, this still leaves us with the choice of which 
three of the kp should be loop momenta. It proves convenient to make different choices 
for different elementary sampling methods. We specify the choice by specifying a triplet of 
integers {(5(1), Q(2), (5(3)} such that / j is fcQ(j). We call Q an index set. Then the complete 
set of propagator matrices can related to the loop momenta by a matrix A: 

kp=j2 Aff^. (18) 

L=l 

Evidently, given the index set Q, the matrix A can be constructed by using the topology of 
the graph. 

Now we characterize certain surfaces, to be called scattering singularity surfaces, near 
which the integrand \f{l) \ is big. To do this, we consider the cuts of our graph in which three 
partons appear in the final state. For each such cut, there is a virtual loop. Let /i be the 
loop momentum. More precisely, li will be the momentum, fcQ(i), of one of the propagators 
in the loop, but we defer for a moment specifying which one. As specified in Sec. |I|, the 
integration contour for li is deformed into the complex plane,Q so that li —>■ /i c = + i^,- 

Before deformation, the integrand is singular on certain surfaces associated with simple 
scattering processes, the scattering singularity surfaces. How can this happen? There are 
four cases. Each case is illustrated by a Feynman graph in Fig. 0. The corresponding 
scattering singularity surface is illustrated in Fig. |[ The cases are 

1. 2 to 2 (s). A virtual parton with momentum li merges with a virtual parton with 
momentum h + h — h to produce a virtual parton with momentum ^2 + ^3- This virtual 



We keep the momenta that enter the final state real. 
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2to2(s) 



2 to 2 (t) 




2 to 3 2to 1 

FIG. 2. Elementary scattering subdiagrams that occur at next-to-leading order. 

— * — * 

parton divides into two partons with momenta I2 and /a that enter the final state. In 
old-fashioned perturbation theory, there is an energy denominator IZ2I + l^sl — |^i| — 
1^2 + ^3 — + which vanishes on the ellipsoid |/i| + 1^2 + ^3 — ^i| = 1^21 + 1^31 in 
the space of li. This is the scattering singularity surface. The contour deformation is 
non-zero on the entire scattering singularity surface. 

2. 2 to 2 (t). A virtual parton with momentum I2 + h scatters from a virtual parton 
with momentum /s — li by exchanging a parton with momentum li. Partons with 
momentum I2 and emerge into the final state. There is a scattering singularity 
surface I/2 + ^i| + 1^3 — ^i| = 1^2 1 + 1^3 1- In this case, there is also a singularity at 
/i = that arises from the propagator of the exchanged parton. This soft exchange 
singularity lies on the scattering singularity surface. Furthermore, in our treatment, 
the contour deformation vanishes at the soft exchange singularity so that, even after 
contour deformation, the integrand is singular there. This is, however, an integrable 
singularity. 

3. 2 to 3. A virtual parton with momentum li collides with a virtual parton with momen- 
tum —li to produce a final state with partons carrying momenta I2, h, and —I2 — h- 
There is an energy denominator IZ2I + l^sl + \h + h\ — 2|/i| + ie, which is singular on 
the sphere = [I/2I + |^3| + 1^2 + ^3]] /2. The contour deformation is non-zero on the 
entire scattering singularity surface. 

4. 2 to 1. A virtual parton with momentum li combines with a virtual parton with 
momentum I2 — h to produce an on shell parton with momentum I2 that enters the 
final state. There is an energy denominator I/2I — — I/2 — li \ +ie, which vanishes on 
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2 to 2 (s) 



2 to 2 (t) 




2 to 1 



FIG. 3. Singularity surfaces associated with the elementary scatterings in Fig. ^. In each case, 
the vectors h and /s (or just h for 2 to 1 scattering) are indicated by arrows. We see the scattering 
singularity surface in the space of h. For 2 to 2 (t) scattering and 2 to 2 (s) scattering, these 
surfaces are ellipsoids. For 2 to 3 scattering, the surface is a sphere, only part of which is shown. 
For 2 to 1 scattering, the surface reduces to a line segment. The integrand is typically not singular 
on the scattering singularity surface because of the contour deformation. However, the contour 
deformation vanishes along the heavy straight lines. Thus, in particular, in the 2 to 2 (t) case the 
integrand is actually singular at li = 0. 



10 



the line li = XI2 with < x < 1. The contour deformation is chosen to vanish at this 
coUinear singularity. As discussed in Ref. @], the collinear singularity cancels when 
one sums over cuts. However, singularities at /i = an li = I2 remain. 

Thus for each cut there are a number of scattering singularity surfaces. There is a contour 
deformation that keeps the integrand from being singular except at special points on these 
surfaces. However, the integrand is sometimes still quite large near these surfaces. For this 
reason, we will design an elementary sampling method for each surface such that the density 
of points is big on the whole surface and singular at the special point if necessary. 

There are two kinds of singularities associated with points kp = where a propaga- 
tor momentum vanishes. First, there is a l/|/cpp singularity when the momentum of the 
exchanged parton m. a, 2 to 2 (t) scattering vanishes. The elementary sampling method 
associated with the 2 to 2 (t) scattering will be designed to take care of this kind of singu- 
larity. Second, there is a l/|fcp| singularity when ant/ propagator momentum kp approaches 
zero. This weak singularity arises simply because massless Lorentz invariant phase space is 
dk/\k\. As it turns out, these singularities in the density will automatically be provided by 
the combined sampling methods without having to specifically arrange for them. 

The 2 to 1 scattering singularity surface is exceptional in the list above in that there 
is no singularity except for the two singular points at /i = and h — h = 0. These two 
singular points are of one or the other of two types mentioned above. Typically a, 2 to 1 
scattering subdiagram is part oi a 2 to 2 (s) 01 2 to 2 (t) scattering subdiagram and the 
two singularities are provided for by the 2 to 2 (s) or 2 to 2 (t) sampling methods. The 
exception is in the case of a self-energy subgraph that is connected to a final state parton. 
In this case, the 2 to 2 (t), 2 to 2 (s), and 2 to 3 sampling methods do not apply and we 
need an explicit 2 to 1 sampling method. Thus we will apply a 2 to 1 sampling method only 
in the case of a self-energy subgraph connected to a final state parton. 

The previous argument indicates that for each scattering singularity surface in the space 
of the loop momentum li in a virtual loop, we should associate a method for choosing li 
that puts a high density of points near this surface. It is then useful to choose the other two 
loop momenta to be the momenta of two of the three partons in the final state. Thus the 
momenta of the final state partons are I2, h, and —I2 — h- The integrand will be singular 
whenever the three partons approach a two jet configuration. Thus we choose the points 
{I2, h} so that the density of points is appropriately singular at the two jet configurations. 

Thus we have a general scheme for organizing the elementary sampling methods. First, 
find all of the possible three parton cuts of the graph in question. Then, for each such final 
state cut, enumerate the scattering singularity surfaces that can occur, counting the 2 to 1 
case only when the virtual loop is a self-energy diagram connected to a final state parton. 
There are five basic situations, as illustrated in Fig. Each combination of a final state cut 
and a scattering singularity surface will correspond to an elementary sampling method. 

One point should be emphasized for clarity. In the numerical method discussed in this 
paper, for each integration point, we compute f/p where / is the integrand and p is the 
density of integration points. Contributions from all cuts of a given graph are calculated and 
summed to form the integrand /. The density p is also a sum, with terms corresponding to 
each of the possible cuts of the graph. Thus there are independent sums over possible cuts 
in / and in p. 
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2tol 2to2(s) 




FIG. 4. Matching of scattering singularities to the structure of one loop virtual subdiagrams. 
A scattering singularity occurs when the energy of an intermediate state matches the energy of 
the final state. For each diagram, the relevant intermediate states are marked with a line through 
the graph. The label near the line indicates the type of the corresponding singularity. For some 
graphs, there is more than one scattering singularity, as indicated. The 2 to 1 singularity is marked 
only in the case of a self-energy subdiagram connected to a final state parton. 
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FIG. 5. Labelling of momenta for graphs of type (d) in Fig. 



Consider the five basic amplitudes with a virtual loop that are illustrated in Fig. ^. In 
two of these cases, there are more than one possible intermediate state involving the virtual 
loop and thus more than one scattering singularity surface. In these cases, it is important 
to understand how the scattering singularity surfaces fit together. 

Let us examine, then, case (d) in Fig. ^ Label the momenta as in Fig. |^. In Fig. ^ we see 
the scattering singularity surfaces in the space of the loop momentum / = li for an arbitrary 
fixed choice of the final state momenta pj. (We draw the figure for / in the plane of the pj). 
The vectors pi, p2, and p^ are indicated. The integration contour is deformed everywhere 
except along the line / = — xpa with < x < 1, which is indicated as a solid line. There is 
an ellipsoidal 2 to 2 (s) surface, |/| + \pi +P2 — 1\ = \pi\ + \p2\- There is also a spherical 3 to 
2 surface, 2|/| = \pi\ + \p2\ + Ipsj. The integration contour is deformed everywhere on both 
of these surfaces, so that the integrand is not singular anywhere. On the other hand, the 
deformation vanishes on the solid line, which can be rather near the ellipsoidal surface in 
the case that the angle between pi and p2 is small. Thus the integrand can be rather big on 
the ellipsoid in this circumstance. Furthermore, the size of the integrand is enhanced where 
the ellipsoid is tangent to the spherical singularity surface. We will need to put an extra 
density of points in the region of this point of tangency. 

Case (e) in Fig. ^ just a little more complicated. Label the momenta as in Fig. 0. Fig. |^ 
shows the scattering singularity surfaces in the space of the loop momentum / for fixed final 
state momenta pj. The vectors pi, p2, and ps are indicated. The integration contour is 
deformed everywhere except along the lines / = xpi, I — pi = xp2 and I — Pi — P2 = xps, 
with < x < 1 in each case. These lines are are indicated as solid lines that form a 
triangle. There are two ellipsoidal 2 to 2 (t) surfaces, |/| + \pi + P2 — l\ = |pi| + \P2\ and 
1^ + 1^1 = IP2I + IpsI- There is also a spherical 3 to 2 surface, 2\l\ = \pi\ + \p2\ + \p3\. As 
in the previous case, (rf), we need an especially high density of integration points near where 
an ellipsoid is tangent to the sphere in the case that this point is near to the triangle, where 
the deformation vanishes. We have also the new feature that each of the ellipsoids shares 
a point with the triangle. At this point the contour deformation is zero, so the integrand 
is actually singular. This is the point where the momentum of the exchanged parton in the 
associated 2 to 2 (t) scattering vanishes. The density of integration points will have to have 
a corresponding singularity at these points. 

This completes the description of the general organization of the sampling scheme. In 
the following sections, we outline its component parts. 
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FIG. 6. Singularity surfaces for graphs of type (d) in Fig. |^. 




FIG. 7. Labelling of momenta for graphs of type (e) in Fig. ^. 



















^ nPi ^^^^ J 





FIG. 8. Singularity surfaces for graphs of type (e) in Fig. ^. 
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V. SAMPLING FOR THE FINAL STATE 



As discussed in the previous section, for each final state cut with three partons in the 
final state, we will define a number of elementary sampling methods. Part of each method 
is to choose with an appropriate density the momentum li that circles through the virtual 
loop that occurs when there are three final state partons. The other part is to choose the 
momenta I2, I3, and —I2 — I3 of the final state partons. We address the choice of the final 
state momenta in this section. 

Let the final state parton momenta be pi, p2, and ps, with J^iPi = 0. (Two of these 
will then be I2 and ^3, but it doesn't matter which ones.) Define -Emax = h^ilPil 

= \Pi\/Ej^ax- Then Z^jX, = 2 and < Xj < 1. We can write the integration over final 
state momenta as follows, using p^ = pi + P2 + 2piP2 cos 9^ so that (icos6'i2 = Pzdp^/ {piP2): 

T = j dpi dp2 dps SC^pi) f{p) 

= J 2El^^^dln E,n^^d cos 9i dcpi d(f)i2 X1X2X3 dxi dx2 f{p). (19) 

Here {9i, 0i} are the angles of pi and 0i2 is the angle between pi x p2 and pi x e, where e 
is an arbitrary reference vector. 

If we sample {Ini^max, cos 6*1, 0i, 0i2} with a density 

_ dN 

d In -Emax d COS 9i d(j)i d(j)i2 

and we sample {xi,X2} with a density 
then the net density is 

p^^^= P^P^ . (22) 
dpi dp2 2E^^^ X1X2X3 

On symmetry grounds, pa should be independent of the angles. Its dependence on -E'max is 
not too important. A convenient choice, normalized so that / dN = 1, is 



1 Stt 



2 



(-E'o/-^max)^ ^ + {Euia.^/ EqY , (23) 



Pa 3 



\3 



where Eq is a fixed parameter with dimensions of mass. 

What should be the properties of ps? We can expect the integrand f{p) to be singular 
in any two jet region: when any of the pj vanishes {xj — > 0) or when two of them become 
collinear (if pi and p2 are nearly coUinear, then X3 ~ 1). The factor l/a;ia;2a;3 in Eq. (p2D 
provides a weak singularity in the Xi ^ regions. We can build in a weak singularity in the 
collinear regions by choosing 

A4 maxj VI - Xj , . 

^ n,vT^ ■ '''' 
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The Xj — i> regions are at the intersections of — 1 regions, so these colhnear singularities 
also enhance the soft singularity of the density. The normalization constant is fixed so that 
/ dN = 1 and is = l/(6-\/2 arcsin(l/-\/3)). One can, of course, change l/i/l — x to 
1/(1 — x)"^ for any A smaller than 1. 

It is a simple matter to actually choose points with the density p in Eq. (p2D. Readers 
interested in the implementation details may consult the notes @] that accompany the 
program code. 

I can comment here about the singularities of the integrand, /. If, for the moment, we 
set the measurement function to 1, then the nominal behavior of / is / oc l/x| for Xj —>■ 
and / oc 1/(1 — Xj) for Xj 1. This nominal behavior is what one gets after performing 
the virtual loop integral, Jdli, inside the {^1,^2} integral |^. These singularities make the 
integration logarithmically divergent. Then the measurement function that is included in / 
is required to vanish in the two jet limit, so that the integration becomes convergent. Now, 
in fact, if we look at the singularities in the {pi,P2} integral before performing the virtual 
loop integral, the singularities are worse than this nominal behavior. For this reason, the 
suppression of the two jet region coming from the measurement function must be suitably 
strong in order to obtain good convergence. It remains for future research to arrange the 
calculation in such a way that the nominal behavior of the integrand / as a function of 
{pi,P2} is obtained. 



VI. SAMPLING FOR 2 TO 2 (S) SCATTERING 



In this section, we consider the sampling method associated with 2 to 2 (s) scattering. 
We need to choose points / = li appropriate to the following case: 

A virtual parton with momentum I merges with a virtual parton with momentum 
^2 + ^3 — ^ to produce a virtual parton with momentum h + h- This virtual parton 
divides into two partons with momenta I2 and that enter the final state. 

In this case, the scattering surface is the ellipsoid \l\ + IZ2 + h ~ l\ = V^l + \h\- The density 
of points should be large on this surface. In order to accomplish this, we use elliptical 
coordinates. 



A. Elliptical coordinates 

The elliptical coordinates ^4+, are defined as follows. First, define k by 

2« = ir2+r3|. (25) 

Now, define coordinates A± by 

A± = ^ (\l\±\L-h-h\). (26) 

Then 
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1 < A+, -1 < A_ < 1. (27) 

The coordinate is constant on elliptical surfaces, while A_ is constant on orthogonal 
hyperbolic surfaces. The scattering singularity is the ellipsoid 

A+ = S+, (28) 

where 

^ ± Ifsl) . (29) 



We need one more coordinate. Let (p be the azimuthal angle of Z in a coordinate system 
in which the z axis lies in the direction of I2 + 13 and I2 has zero y component and a positive 
X component. To state this precisely, define unit vectors 

h + h 



Uy 



h + ^3 
^3 X h 



fix = fly X n^. (30) 

Then 

= arctan (j[ ■ fly /I ■ fi^ . (31) 
The transformation from A^, A_, to / is 

/ = ^(^2 + ^3) + h cos (p fix + It sin (j) fiy + z fi^, (32) 

where 

iT = 4Ai-if^\i-Aiy^' 

z = kA+A_. (33) 
A straightforward calculation shows that the jacobian of the transformation is given by 

dl= dA±^^^ (34) 

PAA(j> 

where 

J— = ,^\Al-Al). (35) 

PAA4, 

The factor {A\ — A^_) is convenient. It provides weak singularities in the density of points 
when A^ ± A_ 0, which corresponds to / ^ and / — /2 — ^3 — ^ 0. 
If we sample points in the variables {v4+, y4_, 0} with a density 

then the total density of points will be 

p = p' X paa^p, (37) 
where p^A^ is given in Eq. (^5]). We next address the question of what to choose for p'. 
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B. The choice of y4 and (f) 



With what density 

^ ^ dA+ dA_ d(t) 

should we sample points in the variables {y4+,y4_,0}? There is no unique answer, but 
consider a density function of the form 

Here r and A are parameters to be specified and the normalization M is given by 



1 27r ^(}/X±l±l 




With this ansatz, there is a high density of points near the scattering singularity surface, 
y4+ = 5*+. The width of the peak is A. Thus A should be matched to the width of the 
peak in the integrand. If the singularity surface is far from the line ^4+ = 1, the contour 
deformation is substantial and peak is broad. On the other hand, if the singularity surface 
is near to the line ^4+ = 1, the contour deformation is small and peak is narrow. I estimate 
that the width of the peak in the integrand is of order (5+ — 1)^ for 5*+ near 1. For large 
5*+ it seems reasonable to choose a width of order 5*+. Thus 1 take 

A = ^-^^j^. (41) 

With the function p' in Eq. (PP|), there is an enhancement of the density of points near 
A_ = ±1. This enhancement is built into the density in order to provide an extra density of 
points near where the ellipsoid in Fig. ^ is tangent to the sphere. As explained in Sec. |V], 
we need a high density of points here when I2 and Z3 are nearly coUinear, that is when 5*+ is 
close to 1. I arrange for this by taking the parameter r in Eq. (|39|) to be 

(42) 

It is easy to sample points {^4+, A., 0} with the density p' in Eq. (BOf) as explained in 



VII. SAMPLING FOR 2 TO 2 (T) SCATTERING 



In this section, we consider the sampling method associated with 2 to 2 (t) scattering. 
We need to choose points / = li appropriate to the following case: 

A virtual parton with momentum I2 + / scatters from a virtual parton with 
momentum ^3 — / by exchanging a parton with momentum /. Partons with 
momentum I2 and /s emerge into the final state. 

In this case, there is a scattering singularity surface 1^2 + /| + I/3 — ^| = \h\ + \h\- There 
is a singularity on this surface at the point where the momentum / of the exchanged parton 
vanishes. The density of points should have a matching singularity at this point. 
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A. Elliptical coordinates 

As in the previous section, we use elliptical coordinates A^,A_,(j). The role of / in the 
previous section is now played by / + /2, so the formulas are a little different. We again define 
K hj 2k = + h\- Now, define coordinates A± by 

^± = ^ {h + l\± \h - /I) . (43) 
Then 1 < and -1 < A_ < 1. Define 

S± = ^ {ikl ± Ifsl) . (44) 
The scattering singularity is the ellipsoid 

= 5+ (45) 

while the soft exchange singularity is at 

{A^,A^} = {S^,S^}. (46) 

We need one more coordinate, an azimuthal angle 0. We define unit vectors {n^,ny,n2} 
according to Eq. (^0|) . Then we define 

(f) = arctan (^{l2 + 1) ■ Uy/ (h + l) ■ ■ (47) 

The jacobian of the transformation from / to A„, 0} is again 1/ pAA<f> where Paa4> is 
given in Eq. (|35|) . If we sample points in the parameters A_, 0} with a density p' then 
the total density of points will be p = p' x pAA<t>- We next address the question of what to 
choose for p' . 



B. The choice of yl+, A_ and 
In this subsection we specify a choice for the density 

^ ~ dA+ dA_ d(j) ^ 

with which we sample in the variables {A^, A_, 0}. We begin by recognizing that we face 
two challenges. First, we would like to have a concentration of points near the surface 
A^ = S+ with an especially high density near A_ = ±1 if 5+ is near 1, as discussed with 
respect to the sampling for 2 to 2 (s). The second challenge is to make p' appropriately 
singular at {A+, 0} = {S+, 0}. We thus take p' to have two parts 

P = a2t Ps + (1 - a2t)pu (49) 
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where a2t is a fixed parameter with < a2t < 1- We take ps to be given by the density (pOD 
for the samphng for 2 to 2 (s). Then we have met the first challenge. It remains to design 
Pt to address the second challenge. 

The main idea is that there is a scattering singularity surface at Aj^ = S^, so that the 
integrand has a factor 

where r) is the amount of deformation of the 5*+ contour. The amount of contour deformation 
vanishes at the soft exchange singularity at {A^, A_,(j)} = {5"+, S'_, 0} and, near to this 
point, 1] is proportional to the square of the distance to {S'+,S'_,0}. Thus, on the surface 
y4+ = S+ near to the soft exchange singularity, we can estimate rj as where 

uj=\A^-S^\ + \(I)\/tt. (51) 

There is another factor of 1/u; that arises from the propagator of the exchanged parton after 
we take real- virtual cancellations into account, as explained in detail in Ref. Thus the 
absolute value of the integrand has a factor that can be estimated as 

for -C 1 and cu^ <^ \A^ — S+\ <^uj. We want pt to have a singularity of the same nature, 
so that the integrand divided by the density of points is singularity free. Thus we take 

A 

^* " [\A+ - S+\+f32tS+uj^]\\A^ - S+\+f32tl2tS+u] ■ ^^^^ 

Here /32t and 724 are fixed parameters and ^ is a rather complicated function of A_ and 
(described below) that is of secondary importance. The main point is that pt behaves like 
(0), with cutoffs when \Aj^ — S+l gets to be smaller than uj"^ or larger than uj. 

We have thus given pt an actual singularity at the point 0} = {5"+, 0}. This 

singularity in pt has a special structure such that matches the structure of the integrand 
/, so that f / Pt is finite in the neighborhood of the soft exchange singularity. This point, 
which is important for the convergence of the numerical integration, was described in some 
detail in Ref. 0. In that paper, however, the construction was not as nice as one would like 
because the ellipsoidal coordinates were not used and the "ridgeline" of p was placed on the 
tangent plane to the ellipsoid rather than on the ellipsoid itself. 

To sample points with the density pt of Eq. (|53|) , we sample first in with a density 
proportional to ln^(7r72f/|0|). Then, with chosen, we sample in A- with a density propor- 
tional to log(72t/ct;)/ct;. Finally, having chosen and A_, we sample in with a density 
proportional to the right hand side of Eq. (|53|). Taking into account the normalization 
conditions, this gives the result in Eq. (|53D with 

^ [ln2(72,) + 21n(72*) + 2] ^ 



A S+P 72t - UJ 



21n(72t/u;) \S+ - 1 + (32tl2t S+uj ^ 



21n^(7r72i/|0|) 



In^ ]+ln^' 



5_ + |0|/7r/ \l-5_ + |0|/ 



TT 



(54) 
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The function A may be a bit complicated, but it is quite benign. In particular, A is not 
singular or zero anywhere provided that > uj holds everywhere. This requires that 
72t > 3. 



VIII. SAMPLING FOR 2 TO 3 SCATTERING 

In this section, we consider the sampling method for the third type of scattering singu- 
larity surface. We need to choose points / = li appropriate to the following case: 

A virtual parton with momentum / collides with a virtual parton with momentum 
— / to produce a final state with partons carrying momenta I2, h, and —I2 — h- 

In this case, the scattering surface is the sphere 2\l\ = I/2I + l^sl + + ^sl- The density of 
points should be large on this surface. There is no special point embedded in the surface 
where the density of points needs to be singular. However, the density should be enhanced 
near the point of this sphere where it intersects the ellipsoid in Fig. |^ or the narrower of the 
two ellipsoids in Fig. ||. 

Since the scattering singularity surface is a sphere, we use spherical polar coordinates 
{r, cos^, 0} to describe /. We choose the 6 = axis along a certain direction Pc- In the case 
(d) illustrated in Fig. ^ we define Pc = —ps- In the case (e) illustrated in Fig. |^, we define 
Pc = Pi if IpiI > IpsI or Pc = -ps if IpsI > \Pi\- 

The jacobian of the transformation from / to {r, cos^, (f)} is 

dr d cos 9 d(f) . . 

dl = , (55) 

Prdij, 



where 



PrH = ^- (56) 



If we sample points in the variables {r, cos^, 0} with a density 

p = -ri — TJi^ (57) 

dr d cos U d(p 

then the total density of points will be 

p = p' X pre4,. (58) 
We now address the question of how to choose p' . Our analysis follows closely that in 



Sec. [VII B| . The main idea is that there is a scattering singularity surface at r = S", where 

5 = l(ir2i + ir3i + ir2+r3i). (59) 

Thus the integrand has a factor 

. ,c^/,^. , (60) 
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where r) is the amount of deformation of the r/S contour. The amount of contour deforma- 
tion vanishes quadratically as {r, cos6'} approaches {Sc, 1}, where 



c 



\P< 



c\ 



(61) 



The point {Sc, 1} is the point / = Pc- That is, / = —p^ in Fig. ^jor in Fig. ^ when |p3| > \pi\. 
Thus, on the surface r = S, the amount of deformation r] can be estimated as where 



COS! 



A = y{i-Sc/sy + {i 

Thus the absolute value of the integrand has a factor that can be estimated as 

1 



\{r/S) 



(62) 



(63) 



for A ^ 1 and A^ <^ \{r / S) — 1\ ^ A. We want p' to have a singularity of the same nature, 
so that the integrand divided by the density of points is singularity free. Furthermore, 
we would like p' to have an extra factor of 1/A^ to cancel the factor from the singularity 
associated with the narrow ellipsoid in Fig. ^ as we approach ^ = 0. Thus we take 



B 



P 



A (|(r/5) - 1| + a3A2) {\{r/S) - 1| + a^A) ' 



(64) 



where as is a fixed parameter and where ;B is a rather complicated but nonsingular function 
that gets the normalization right: 



B 



asjl-A) 
ISttS 



/ 1 l + asA^ \ 
""[a^ 1 + asA J 



In 1 + 



:i-sc/s)\ 



(65) 



IX. SAMPLING FOR 2 TO 1 SCATTERING 

In this section, we consider the sampling method for the fourth type of scattering singu- 
larity surface: 

A virtual parton with momentum / combines with a virtual parton with momen- 
tum ^2 — ^ to produce an on shell parton with momentum I2 that enters the final 
state. 



As mentioned in Sec. I^, the 2 to 1 scattering singularity surface is exceptional in that there 
is no singularity except for the two singular points at / = and /2 — ^ = 0. Typically a 2 
to 1 scattering subdiagram is part oi a. 2 to 2 (s) oi 2 to 2 (t) scattering subdiagram and 
the two singularities are provided for by the 2 to 2 (s) or 2 to 2 (t) sampling methods. The 
exception is in the case of a self-energy subgraph that is connected to a final state parton. 
In this case, the 2 to 2 (t), 2 to 2 (s), and 2 to 3 sampling methods do not apply and we 
need an explicit 2 to 1 sampling method. Thus we apply a 2 to 1 sampling method only in 
the case of a self-energy subgraph connected to a final state parton. 
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We will choose / with a density that is a linear combination of five densities: 

P(0 = «laPa(0 + (^lbPb{l) + OiiaPa{l " ^2) + Oilb Pb{l - h) + (1 " 2aia - 2aib)pc{l)- (66) 

Here aia and aih are fixed positive parameters with (1 — 2aia — 2aib) also positive. 
For Pail) we take 

The density ph is a simple variation on this with IZ2I replaced by a scale M: 

4^ l^{\l\^M\ 



Here M = (I/2I + l^sl + K2 + ^sD/S, where the momenta of the particles in the final state are 
li, I3, and —I2 — h- The pa and pb terms provide singularities at P = and {I2 — l)"^ = 0. 

The pc term provides a non-singular density of points in the neighborhood of the collinear 
line / = x/2 with < x < 1. For Pc(0 take 

Pcil) = , , ^ 2- (69) 

2^ v//3? + (x- 1/2)2 (^^2 + (^_i/2)2 + (fi + ^2^;2) 

Here /?i and 71 are fixed parameters and x and It are defined by writing 

1 = xh + fr, (70) 

— * — # — * 

where It is the part of / that is orthogonal to l2- 

X. CONCLUSIONS 

In a next-to-leading order calculation of three jet cross sections and similar observables 
in electron-positron annihilation, a given graph leads to several contributions to the cross 
section. Each of these contributions has the form of a measurement function times a quantum 
amplitude times a complex conjugate amplitude, as in Fig. We must integrate the sum 
of these contributions over three loop momenta. In this paper, I have presented a method 
for sampling the space of loop momenta in order to perform the integrations by the Monte 
Carlo method. 

The main organizing principle for this sampling is to construct the density of integration 
points as a sum and to have one or more terms in this sum correspond to each possible three 
parton final state for the graph. 

A cut graph with a three parton final state has a virtual loop. We are thus led to consider 
a simple problem in quantum mechanics: where in the space of the loop momentum are the 
singularities for two body to n body scatterings? The generic answer is that the singularities 
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lie on ellipsoidal surfaces. This leads us to choose elliptical/hyperbolic coordinates for the 
loop momentum space. For the most part, we avoid letting the integrand be truly singular 
at these surfaces by giving the loop momentum a suitable imaginary part according to the 
contour deformation recipe of Ref. [Q. Nevertheless, it is helpful to make the density of 
integration points large near these ellipsoidal surfaces. 

We are left with special points on the scattering singularity surfaces where the contour 
deformation vanishes, so that the integrand is actually singular. This singularity corresponds 
to soft parton exchange and is significant in gauge field theories. We have seen how to give 
the density of points a singularity structure that matches that of the integrand as one 
approaches a soft singularity. 

The sampling method that has been described here is far from optimal. It is, however, 
at least reasonably systematic, and it gives good results for three jet observables that do not 
get significant contributions from parton final states that are near to two jet configurations. 
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